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ABSTRACT 


This publication presents a novel multistage estimation scheme for estimating the 
parameters of a received carrier signal possibly phase-modulated by unknown data, and 
experiencing very high Doppler, Doppler rate, etc. Such a situation arises, for example, 
in the case of Global Positioning Systems (GPS) where the signal parameters are directly 
related to the position, velocity, acceleration and jerk of the GPS receiver. 

In the proposed multistage scheme, the first stage estimator operates as a coarse 
estimator resulting in higher rms estimation errors but with a relatively small proba- 
bility of the frequency estimation error exceeding one-half of the sampling frequency 
(an event termed cycle slip). The second stage of the estimator operates on the error 
signal available from the first stage, refining the overall estimates, and in the process 
also reduces the number of cycle slips. The first stage algorithm is selected to be a 
modified least squares algorithm operating upon the differential signal model and re- 
ferred to as differential least squares (DLS). This estimation stage provides relatively 
coarse estimates of the frequency and its derivatives. The second algorithm is simply 
an extended Kalman filter (EKF) which also yields the estimate of the phase along with 
a more refined estimate of frequency as well. 

A major advantage of the proposed algorithm is a reduction in the threshold on 
received carrier power-to-noise power spectral density ratio (CNR) as compared to the 
threshold achievable by either of these algorithms alone. In fact, it appears from the 
simulations that for the case of an unmodulated carrier, the proposed scheme achieves 
the same threshold as for an almost exact and computationally intensive implementation 
of the maximum likelihood estimator (MLE). For the case when data modulation is 
present, the proposed scheme provides an improvement of about 6 dB in terms of CNR 
compared to an earlier MLE scheme reported in the literature. The overall complexity 
of the algorithm is about two times the complexity of a third-order Kalman filter or a 
single fourth-order EKF. 
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1. INTRODUCTION 


The problem of estimating the parameters of a received quasi-sinusoidal signal in 
the presence of noise occurs in diverse scientific and engineering disciplines [1-15]. The 
signal parameters of interest are usually the phase, frequency and frequency deriva- 
tives which are varying with time. The estimation problem becomes considerably more 
difficult if the received carrier is modulated by unknown data while simultaneously 
experiencing very high Doppler and Doppler rate. Such situations occur, for exam- 
ple, in the cases of Global Positioning System (GPS) receivers and NASA deep space 
communication links under high spacecraft dynamics. 

In a previous publication [7], an estimator structure based on the maximum like- 
lihood estimation (MLE) of code delay and Doppler frequency over a single data-bit 
period has been proposed and analyzed for the GPS applications. This scheme esti- 
mates Doppler frequency (assumed constant) over successive intervals of bit periods, 
followed by a Kalman filter tracking Doppler frequency and frequency rate. The scheme 
does not involve the carrier phase estimation. For the dynamic trajectories simulated in 
[7], the approximate ML estimator performance exhibited a threshold of about 30 dB- 
Hz in terms of the received carrier power-to-noise power spectral density ratio ( P/N 0 ), 
below which rapid performance deterioration occurred. 

For GPS applications, an alternative scheme has been proposed [8] wherein a par- 
allel (non-dynamic) link is established between the GPS satellites and a control ground 
receiver for the purpose of communicating the data to the ground receiver. The ground 
receiver simultaneously receives the frequency- translated version of the GPS receiver 
signal and removes the data modulation from this dynamic signal. Such an effectively 
demodulated signal is then processed by the estimation algorithm to obtain the desired 
signal parameter estimates. There are several estimation schemes in the literature for 
this problem. See for example [9-14]. 

More recently in [9] a scheme for simultaneous detection and estimation has been 
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proposed. This scheme is based upon first estimating the received signal's local (data 
dependent) parameters over two consecutive bit periods, followed by the detection of 
a possible jump in these parameters. The presence of the detected jump signifies a data 
transition which is then removed from the received signal. This effectively demodu- 
lated signal is then processed to provide the estimates of the global (data independent) 
parameters of the signal related to the position, velocity, etc., of the receiver. A key 
feature of this scheme is that to a certain extent the data detection is independent of 
the acquisition of the phase or frequency of the received carrier signal in contrast to the 
conventional decision-directed phase-locked loop receivers. In these latter receivers, the 
data detector is an integral part of the loop and depends upon the acquisition of the car- 
rier phase and/or frequency. Thus, under low CNR and/or high dynamic conditions, 
the loop may not acquire lock or frequently lose it during tracking. From the simula- 
tions of [9], it is seen that the scheme offers very significant improvement in terms of 
the required CNR over the AMLE algorithm of [7]. The detection-estimation scheme of 
[9] has a computational complexity about three times that of a single extended Kalman 
filter. 

In this publication we propose an alternative scheme for the estimation of the signal 
parameters applicable for both the cases of unmodulated carrier and when the signal is 
phase-modulated by unknown data. The proposed scheme is somewhat simpler than 
that of [9] in that it is not essential to detect the data modulation explicitly. Basically, the 
proposed algorithm involves an appropriate modification of the DLS scheme of [11] so 
as to apply the algorithm to the case of unknown data modulation. As discussed in [11], 
if the DLS technique is applied with the Nyquist sampling of the received signal, there 
is expected a loss in performance compared to the optimum achievable performance. 
The techniques of [11] propose oversampling and cyclic sampling strategies to avoid 
such a loss. In the present example of GPS application, we sample at Nyquist rate and 
propose an alternative method to keep the overall performance close to optimum. The 
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technique proposed here consists of a multistage procedure wherein the parameters of 
the signal are estimated in more than one stage. First, the parameters are estimated by 
an algorithm like DLS which has a low threshold on CNR but with possibly higher rms 
estimation errors. Then an error signal whose parameters are equal to the difference 
between the true parameters and the above estimates is processed by another algorithm 
to estimate these error signal parameters. Since the error signal involves much smaller 
dynamics, the second algorithm can have smaller bandwidth resulting in a smaller 
estimation error. In principle, this procedure can be repeated any number of times with 
successive algorithms having progressively lower bandwidths. 

In this publication we confine ourselves to two stages of recursion and apply a 
third-order extended Kalman filter algorithm for the second recursion. It is expected 
from such an estimation structure that the overall algorithm would have both smaller 
threshold and a smaller estimation error compared to either algorithm operating by 
itself. Indeed, this is borne out by simulations presented in the publication. Thus, for 
the case of no data modulation, whereas the threshold on SNR is about 1.5 dB lower than 
the extended Kalman filter, the estimation errors are only marginally higher than for the 
third-order EKF alone. The threshold achieved is in fact the same as that achieved for 
a nearly exact implementation of the maximum-likelihood estimator (MLE). It is also 
noted that the threshold achieved by the proposed scheme is about 3 dB lower than 
conventional cross product AFC (CPAFC) loops and phase-locked loops [14], whereas 
the rms error is less than one-half of that obtained by CPAFC. The rms error is marginally 
higher than EKF due to the non-optimal sampling used in the DLS algorithm. 

For the case of data modulation, we compare our results with those of [7], where 
analysis and simulations are presented on the performance of Fast Fourier Transform 
(FFT) based MLE algorithm. In [7], the trajectories of the GPS signals have somewhat 
less severe dynamics compared to those considered in this publication. In terms of 
threshold on CNR, the proposed scheme of this publication exhibits a threshold of 24 
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dB-Hz compared to about 30 dB-Hz reported in [7], thus providing an improvement of 
about 6 dB. In terms of the rms frequency estimation errors, at a 30 dB-Hz CNR, the 
scheme of [7] provides a rms range rate error of about 6 m/s compared to an error of less 
than 2 m/s achieved in this publication. There is also a very significant improvement in 
terms of the rms position estimation error. At about 30 dB-Hz an rms error of 1 meter 
is reported in [7], compared to about 0.25 meter obtained by the proposed algorithm. 
It may also be remarked that in the previous scheme, pseudo-random codes with rate 
10.23 MHz are needed for the purpose of range measurements, thus requiring a zero- 
crossing channel bandwidth in excess of 20 MHz. The proposed scheme, on the other 
hand, extracts the range information from the carrier signal itself and thus needs a 
bandwidth equal to only a fraction of 1 MHz. 
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2. THE RECEIVER CONFIGURATION 


We consider the problem of estimating the high dynamic phase process ©/(<) of the 
desired signal s/(t) observed in the presence of an additive narrow-band noise process 
vj(t) as 


ri(t) = s/(t) + v/(t ) 

= ASin(w c t + 0/(t) 4- ir D(t)) + r>/(<) 


( 1 ) 


In (1) is the received signal carrier frequency in the absence of any dynamics and 
D(t) is a binary digital waveform taking on the possible values of 0 or 1. In the case 
of a Global Positioning System (GPS) receiver, the process ©/(<) arises from the receiver 
dynamics, and over a sufficiently small estimation period. 


©z(t) = 0/o + w/oP -I- jiot 2 + 0/o< 3 (2) 

for some unknown parameter vector V’/o = [0/o w /0 j J0 0/ O ]'. In a somewhat simpler version 
of the problem, the data modulation D(t) is either absent or is assumed known and thus 
can be eliminated from (1). In the sequel both of these cases would be treated in some 
detail. 

In the first instance we only estimate the parameters related to the frequency and 
its derivatives using the DLS algorithm. For this purpose, the received signal r/(t) is 
quadrature demodulated by the voltage-controlled oscillator (VCO) signal s L (t) as shown 
in Figure 1. We assume that the input to the VCO is a signal which is an appropriate 
quadratic function of time t resulting in its instantaneous phase 0 L (t) in (3) below. 


«i(<) = 2 Cos(w c t + ©/,(<)) 

= @L0 + WLO^ + "/Lot 2 + &Lot 3 


(3) 


for some constant vector y> L0 = [0 LO w L0 j L0 <5 L0 ]'. The sampled version of the in-phase 
and quadrature components of the demodulated signal are given by 
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y(k ) = ASin(Q(k) + wD(k)) + t>,(fc) 

z(k) = ACos(Q(k) + irD(k)) + v q (k) ; k = 1, 2, . . . , N 

where 


( 4 ) 


0(*) = 0/(k) - e L (k) = 6 0 +Lj 0 kT,+ 70 (kT,) 2 +6 0 (kT.) 3 
V’O — V’/O — V’LO = [*0 Wo 7o <5o] / 

and V'o is the parameter vector characterizing the error signal to be estimated, with T, 
denoting the sampling interval. In (4) «,•(*) and v q (k) represent the sampled in-phase and 
quadrature components of the bandpass noise process v I (k). The parameter vector V'o is 
estimated by the Differential Least Squares (DLS) algorithm of [11], as described in the 
following section. 
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3. DIFFERENTIAL LEAST SQUARES (DLS) ALGORITHM 

Consider first the problem of estimating the unknown parameters w 0 , y 0 and 6 0 from 
the measurement (4) for the case of no data modulation, i.e. when D(k) = o. Following 
[11] we expand Sin(e(t)) in a Taylor series around **_ i = (k - 1 )T, to obtain 


Sin(0(t)) = Sin{Q(k - 1)) + TQ(k - 1 )Cos(G(k - 1)) + . . . (5) 

with a similar expansion for Cos(0(f)). For small (t - t k ~ x), the series in (5) may be 
approximated by the first two terms and from (4) one obtains the following differential 
signal model with r k = (k- i)T,: 


y d (k) = y(Jfc) - y(k - 1) = T,(u> 0 + 2 To r t + 3 6 0 ^)z(k - 1) + &(Jb) 
z d (k) = z(k) - z(k - 1) = —T, (w 0 + 2 T or t + Z6 0 r^)y(k - 1 )+£,(*) 

where 


( 6 ) 


£i(k) = Vi(k) - Vi(k - 1 ) - T,(u> o + 2y 0 r k + 3 S 0 ^)v t (k - 1 ) 


* f (*) = v t (k) - v q (k - 1) + T,(u> o + 2 7o r k + 35 0 r t >,(fc - 1) 


(7) 


The measurement model (6) may be rewritten in the following standard form as in [11]: 


Z d {k) = H'{k)p + m ;k = 1,2,. ..N 


( 8 ) 


where 


0' = [w 0 270 6(5 0 ] 



T,z{k - 1) 

-T,y(k~ 1) 


T.r k z(k - 1) 
-T,r k y(k-1) 


0.5 T,rj*z(k — 1) ‘ 
—0.5 T,r£y(k — 1) 


(9) 


Z' d {k) = [y d (k) z d (k)] ;£'(*) = Uk) i q (k)] 

The parameter vector 0 in (8) is now estimated by an exponential data-weighted least 


squares algorithm in a recursive or nonrecursive form (Kalman filter). In its nonrecursive 
form the estimate of 0 obtained on the basis of N measurements and denoted by 0{N) 
is given by 
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( 10 ) 


0(N) = 

i = i ;'=i 

where A is some appropriate weighting coefficient with 0 < A < l. An equivalent recursive 
form of (10) is the following algorithm: 

0(k) = p{k - 1) + L(k)e(k) 

L(k) = P(k - l)H(k)[XI + H'(k)P{k - l)tf (fc)]" 1 

( 11 ) 

P(k) = {P(k - 1) - P(k - l)H(k)[XI + H'(k)P(k - l)H(k)]-'H'(k)P(k - 1)}/A 
f (k) = Z d (k) - H'{k)ji{k - 1) ;k = 1,2, ... ,1V 

Note that the matrix to be inverted in (11) is only a (2 x 2) matrix. In an alternative 
but equivalent form one may process the scalar measurements j i d (k), z d (k ) sequentially 
instead of working with vector measurement z d (k). Moreover, the matrix P(k) of (11) 
with k = N is the same as the matrix inverse in (10), i.e., 


k 


3 = 1 


Alternatively, the matrix P~'(k) may be written as 


( 12 ) 


P~\k) = £A ^{z^ + y^BjThBj i 

3 = 1 


?1 


T i 


TJ Tj 

r? r? 

3 3 

r? r 4 
3 3 


(13) 


and thus the matrix P(k)H{k) required in the update of (3(k), and equal to l(Jfc) in (11) 
may be approximated by 


P(k)H(k) S *(*)[*(*) - y(k)] 

X'(k ) = (£ X k ^Bj)-\A 2 + r k rD'Tr 1 

3=1 

where 


( 14 ) 


E[z\j) + y*(j)] = A 2 + al 


®[«ii = = r 
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In (14) the vector x(*) is data independent and thus could be precomputed for k = 
1 , 2 ,..., N for computational simplification of (11). Similarly, in the implementation of 

N 

(10), the first matrix may be replaced by the data independent matrix (£ A N - j B j )(A 2 + 

; = 1 

°l)Tl 

Modified Least Squares Algorithm 

If the noise £(Jfc) in the signal model (8) were white, then the estimate 0(N) obtained 
from the algorithm (10) or (11) would approach 0 as N — oo, if one ignores the approx- 
imation made in arriving at model (8) and A is selected equal to 1. However, as the 
noise £(Jb) in the model (8) is colored, there would be considerable bias in the parameter 
estimates under low to medium signal-to-noise ratios. To reduce such a bias or pos- 
sibly eliminate it we propose the following simple modification. If the instantaneous 
frequency u(r k ) at time r k given by (w 0 + 27 0 r fc + 36 0 r t 2 ) appearing in ( 7 ) is small compared 
to l/T„ then the noise vector £(*) is equal to v(/fc) - v(k - l) where v(k) = [v { (*)v ( (Jr)]'. This 
situation may be depicted by Figure 2a. 

To eliminate the bias, the noise f(fc) must be whitened by passing through the trans- 
fer function (l-z -1 ) -1 as shown in Figure 2b. The least squares algorithm is, in general, 
nonlinear and time-varying. However, if we assume that the algorithm in (11) asymp- 
totically approaches a time-invariant system, then under such an assumption, one may 
interchange the least squares algorithm with the transfer function (l - * -1 ) -1 to arrive 
at Figure 2c. This, of course, corresponds to post-averaging the least squares estimates. 
Such a simple procedure provides very significant improvement in the estimates when 
the signal-to-noise ratio is low. In the simulations of the next section, the infinite time 
averaging is replaced by an exponentially data-weighted averaging to take into account 
the time variation of the parameters to be estimated. In Figures 2b and 2c, p uB denotes 
an unbiased estimate of /?. 
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DLS Algorithm in the Data Modulation Case 

In this case the data samples D(k) in the signal model (4) may take possible values 
±1 and the received signal may equivalently be written as 


y{k) = D(k)ASin(Q(k )) + v,(fc) 


( 15 ) 


z(k) = D(k)ACos(Q(k)) + v g (k) ; k = 1, 2, . . . , N 

Thus, as may be easily verified, over any bit interval 71, where D(k) remains constant, 
the differential model (6) remains valid irrespective of the value of D(Jfe). The model (6), 
however, is not applicable for those samples which lie on the bit boundaries, i.e., when 
y(k) and y(k - l) lie in different bit intervals. A simple modification of the algorithm to 
take care of the data modulation case is to simply discard such differential samples. If 
the number of samples M over any bit period is fairly large, this would incur a negligible 
loss in the effective signal-to-noise ratio compared to the case of no data modulation. 
In fact, such a loss is simply equal to iO/o^ 10 (i - -fa) dB which is 0.45 dB for M - 10. This 
is corroborated by the simulations of the next section. 


Estimation of Time- Varying Parameters 

In the signal model considered above we have assumed that the input signal pa- 
rameter vector rl> 10 is either a constant or a slowly varying function of time. However, 
in practice this may be the case only over relatively short intervals of time, but there 
may be large variations in V/o over a comparatively large observation period. To take 
into account such a variation and to ensure that the instantaneous difference frequency 
n(t) = ^0(7) (the sampled version of 0(<) given in (4)) remains within the low-pass filter- 
pass band of Figure 1, the parameter vector 4> LQ generating the instantaneous frequency 
of the VCO is updated at regular intervals of T = NT, sec for some integer N. The pa- 
rameter vectors i>i, V>lo and rp 0 would change their values at intervals of T sec, assuming 
that the value of N is selected to be sufficiently small so that the variation in ip f over 
any T sec interval is small. Denoting by 0 L o(T+),uj lo (T+), etc., the values of reference 
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oscillator parameters just after the update at time T, we have that 


0lo(T+) = Qlo(T-) 

ulo(T+) = n L0 (T-) + u>o(0/T) + 2Tyo(0/T) + ZT 2 6 o (0/T) 

( 16 ) 

Jlo(T+) = r t0 (T-) + ZT6 o (0/T) 

6 i0 (T+) = A L0 (T~) + 6 o (0/T) 

In equation (16), G L o(t-), ^lo(T-), etc., represent the oscillator instantaneous phase, 
frequency, etc., just before the correction and the remaining terms on the right-hand 
side of (16) represent the correction made on the basis of the estimation algorithm. 
Thus, 

©to (T-) = e L0 ( o+) + w LO (o+)r + 7 £o(o+)r 2 + * io (o+)r 3 

f2io(T— ) = w/;o(0+) + 27 i,o( 0+)T + 35io(0+)T 2 

(17) 

r(T-) = 7 Lo( 0 +) + 3TM0+) 

A L o(T-) = ^ o ( 0 +) 

Note that in (16) there is no correction in the oscillator phase as the DLS algorithm 
does not provide the phase estimate. In equation (16), w 0 (o + /T) denotes the estimate 
of parameter w 0 (o+) obtained on the basis of measurements up to time T. Since there 
is no step change in the oscillator phase, the sampled measurements y(N) and z(n) at 
the demodulator output are the same with or without a correction at the instance NT,. 
However, the subsequent measurements y(N + j) and z(N + j) are now expressed with 
respect to the new parameter vector VolT+j = [0 O (T+) w 0 (T+) 70 (T+) <5 0 (T+)]' = V/o(T+) - 
V , lo(T+) as in (18) below. 

y(N + j) = ASin(Q(N + j)) + v { (N + j) 

z(N + j) = ACos(G(N + j)) + v q (N + j) ; j = (18) 

Q(N + j ) = 0 O (T+) + u>o (T+)jT, + 7o(T+)OT,) 2 + 6 0 (T+)(jT,f 
The last three elements of the vector ip 0 (T+) will be zero if there is no change in the 

input signal parameters over the T sec interval and the estimate of V>o(0+) is obtained 
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with zero estimation error. We thus set the a-priori estimate of the vector 4> 0 (T+) equal 
to 0 and apply the DLS algorithm to estimate Vo (T+) on the basis of observations {y(N + 
j),z(N + j)\j = l, . . . , N}. The measurement model is obtained by simply replacing the 
index k by k + N in y(k),z(k),^(k),^(k) in equations (6-9) but with r k - (k- \)T, as before 
(corresponding to a shift in time reference). 

In the estimation of Vo (T+) via the recursive algorithm (11) with the index k = N + 
l, . . . ,2 N, the "initial" covariance matrix P(N + l) is obtained as 


P(N + 1) = A FP(N)F' + Q 


(19) 


where 


F = 


1 

0 

0 


2 T 3 T 2 
1 3 T 

0 1 


and the matrix Q represents the uncertainty introduced due to the change in the input 
process parameters over the interval of T sec. Specifically, the last diagonal element of Q 
represents the variance of the change in the parameter 6 6 (equal to the second derivative 
of frequency and related to the jerk of the physical trajectory) over the interval T. The 
above procedure is then extended in a straightforward manner to the subsequent update 
intervals. The estimates of input signal phase and frequency at time instances er+ are 
then simply given by 0 LQ {tT-\-) and u> L0 ((T+) respectively for ^ = o, 1,2,. ... 
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4. MULTISTAGE ESTIMATION 


We note that most of the phase and frequency estimation schemes can be repre- 
sented as in Figure 3a. The block representing the VCO correction signal could be 
any recursive or semi-recursive algorithm including an EKF or DLS algorithm, and the 
VCO update interval may be some integer multiple of sampling period T,. Figure 3b 
depicts an equivalent and a more compact representation of the estimation scheme. An 
important observation made from figure 3b is that along with the estimate of the in- 
put phase process © 7 (t) represented by 0 7 (<), there is also available a pair of signals, 
(ASinQ(t), ACosO(t)) dependent upon the estimation error Q(t) = Qi{t)-Q L {t). These error 
signals have exactly the same form as the signals at the input to the estimator. Moreover, 
the additive noise associated with these signals has statistics identical with the statistics 
of the noise at the input to the estimator. Therefore, this leads to the clear possibility of 
estimating the error signal 0(t) in a way similar to the estimation of 0 7 (<). In fact, the 
procedure can at least in principle be repeated any number of times as shown in Figure 
4. It may also be noted that the first estimation stage (Figure 3a) requires a VCO for 
down conversion as the actual input to this stage is at rf frequency u e . However, sub- 
sequent estimator stages generate the error signals by simple baseband computations. 
For example, in the discrete-time version of the estimation procedure, the signal at the 
output of estimation stage m may simply be computed as 


y m (k) = y m - 1 (^)Cos(0 m - 1 (/b - 1)) - z m - l (k)Sin(& m - 1 (k - 1)) 


( 20 ) 

z m (k) = y m ~ 1 (k)Sin(Q m ~ 1 (k - 1)) + * m - 1 (fc)Cos(0 m - 1 (fc - 1)) ; m = 2, 3, . . . n 

Note that 0(t),€,(t) and v q (t) of Figure 3b are respectively equal to 0 HO/ »<(*)/ »J(0 of 
Figure 4. The refined estimate of 0 7 (ifc) is then simply given by 


&i(k) = G L (k) + e l (k) + -. • + ©"(*) (21) 

The advantage of such a recursive estimation procedure is that the overall threshold 
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in terms of the required CNR for the multistage estimator can be made much smaller 
than a single stage estimator, especially in situations involving high dynamics. In a 
single stage estimator, due to the high dynamics involved, the process parameters may 
be assumed to remain constant only over short intervals of time. Thus, the estimator is 
forced to use a relatively large noise bandwidth (shorter averaging period), resulting in 
large errors in the phase and/ or frequency estimates. If the estimation errors are outside 
the region over which the error model (linear) assumed for the estimator remains valid, 
the estimator is said to be working below threshold or in the out-of-lock condition. In 
this condition the estimation errors can be several orders of magnitude higher compared 
to the operation above threshold. In a multistage estimator environment this difficulty 
can be circumvented by successive reduction of the dynamics (the estimation errors 
due to dynamics) at the output of consecutive estimator stages and by averaging the 
signal over progressively longer intervals (and thus progressively reducing the effect 
of noise) over which the process parameters remain nearly constant. In this estimation 
structure, none of the individual stages (except the last one) need necessarily operate 
above its threshold. For the convergence of the overall estimator, one only requires 
that the estimates are made in the right direction (estimation errors do not exceed the 
parameters to be estimated in some average sense). The proposed estimation scheme 
may look familiar if one compares it with the standard practice of weighing wherein 
the estimation proceeds from a coarse estimate to a successively refined one thereby 
achieving high estimation accuracy with a comparatively less complex setup. Similar 
approaches, variously termed method of successive approximation, method of sieves, 
and so on, are also common in various mathematical disciplines including theory of 
differential equations, probability theory, etc. 
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5. RECURSIVE DLS-EKF ALGORITHM 


In this publication we consider a simple special case of the more general multistage 
estimator of Figure 4, wherein the first stage is the DLS algorithm of the previous section 
and the second stage is an extended Kalman filter (EKF). As the dynamic variation of 
the error signal e(k) = Q^k) of Figure 4 is much smaller compared to the original signal 
Q(k) (the frequency variation over any update interval is much smaller), the effective 
averaging time period for the second-stage Kalman filter can be selected to be higher 
than for the first stage DLS algorithm. This is achieved by selecting a smaller value 
of the "dynamic noise" covariance matrix Q (see equation (19)) for the EKF and/or a 
higher value for the exponential data weighting coefficient A. 

Extended Kalman Filter (EKF) Algorithm 

Here we consider the problem of estimating the unknown error signal parameters 
w 0 1 , 70 1 and Sot in the £th VCO update period for any integer l > l on the basis of the set 
of measurements {^(Ar),* 1 ^)} of equation (22) below (see also Figure 4) by an EKF. 


y 1 (k) = ASin(Q 1 (k)) + t >}(k) 
z 1 (k) = ACos(Q\k )) + Vq (Ar) 

= 6 oi + uotjT, + yoi(jT)) 2 + 6ot(jT,) 3 


( 22 ) 


k = N(e-l) + j-,j=l,2,...,N ; £= 1,2, . . . 

Note that as for the DLS algorithm the parameter vector i> ol = [0 O/ w 0/ y ot <W]' may 
be different over different VCO update intervals. For computational simplicity we use 
a third-order EKF and thus ignore the contribution of the last term in the expression 
for ©Hfc) which is appropriate for the GPS trajectories considered here. Denoting the 
state and parameter vectors at time k = N(£- l) + j by Mi) and % respectively, i.e., with 
ip(j) = [l jT, 0.50'T,) 2 ]', T) t = [dot w 0 1 270 /]', the extended Kalman filter equations for 
the update of %, the estimate of are given by 
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mU) = m(j - 1) + Kt{j)v t {j) 

Ki{i) = E,(i - !)iKj)( A + ^'U) E/ J ’ ~ 1 W))” 1 

!>•) = <E 4 o- - 1) - E/j - wow + v-'O) E<o - war vco E<o - m/* ( 23 ) 

Vt(j) = y^^Cos^O - )) - ^(^O^OM.?)) 

Mi) = tf'OJfeO - 1) ; * = (e - 1)N + j ; j = 1 , 2 . . . , N; t = 1, 2, . . . 

In the equations (23) above, the initial estimate rj t { 0) is simply taken to be equal to 

fn-i(N). This is an appropriate choice for the initial estimate in view of the fact that if 

the first stage of the estimation algorithm is convergent then Tj t would be random and 

remain close to zero for all L However, if the errors in the previous stage are not close 

to zero, then the components of % would possess some continuous drift term, i.e., u 0i 

will have a component linear in time if jot ^ 0. The "initial error covariance" matrix 

X^(o) is simply set equal to some diagonal matrix representing the uncertainty in the 

difference parameter ru -rjt-i. 

Estimation in the Presence of Data Modulation 

In this case one could apply the more sophisticated version of [9], wherein an explicit 
detection of possible data transitions is followed by the demodulation of data, thus 
effectively reducing the problem to the case of no data modulation considered above. 
However, here we bypass such a detection and instead propose a simple modification 
in the estimation algorithm that takes into account the data modulation. If the VCO 
update interval T is selected equal to bit period T b , then the data modulation represents 
an additional phase uncertainty at the boundaries of the update intervals. This is taken 
into account by adding an appropriate value, say (?r/2) 2 , to the first diagonal element of 
y^(0) and modifying the initial estimate %(0) by tt/ 2, i.e., %( 0) = fu-i(N) + tt/ 2, for those 
values of t that correspond to bit boundaries. Such an algorithm is expected to result in 
somewhat higher estimation errors compared to the more sophisticated scheme of [9], 
but is much simpler in terms of implementation. In this case of a two-stage estimator 
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(see Figure 4), the estimates of the input signal phase and frequency at time instances 
IT are given by 


4>(£T) = e L ((t - 1 )T+) + 4 >'(n)th(n) 

= Q t ((t - l)T+) + u, ol + 2w • (NT,) 

where w ol and represent the second and third element respectively of m(N). 


( 24 ) 
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6. SIMULATIONS 


In the following we present simulation results obtained when the algorithm is ap- 
plied to the tracking of phase and frequency for high dynamic GPS receivers [14]. For 
the purposes of simulation we assume that the pseudo-random code has been removed 
from the received signal and symbol timing has been acquired. We consider both the 
cases when the data modulation has also been removed via an auxiliary link [8, 14] 
and the case when an unknown data modulation is present. For the simulations we 
consider a sampling rate of 500 samples /second and simulate a high dynamic trajectory 
considered previously in [14] and reproduced in Figure 5. In Figure 5, the acceleration 
and the jerk (the derivative of acceleration) are measured in units of g where g is the 
gravitational constant (equal to 9.8 m/s). In the case when data modulation is present, 
we assume a BPSK modulation at a rate of 50 bits/second. 

The parameters of most interest in this application are the instantaneous phase and 
frequency of the input signal rj(t), which corresponds to the high dynamic GPS trajectory 
of Figure 5. Since we are mainly interested in the tracking performance of the proposed 
algorithm, we assume as in [14] that the initial trajectory parameters at zero time are 
known. The received signal carrier frequency f c = w c / 27 r in the signal model (1) is taken 
to be equal to 1.575 GHz. The GPS receiver instantaneous pseudo-range R in meters and 
velocity v d in m/s are related to the instantaneous phase 0 7 (/) of (1) and its derivative 
0/(f) as 


/?= — — 

2* fc 

Vd = fdj- = 


e_c_ 
2tt f c 


( 25 ) 


where f d denotes the Doppler in Hz and c is the speed of light. Denoting by R(() and 
f d {l) the estimates for R(t) and f d (£) respectively, which denote the range and Doppler of 
the input trajectory at the end of «h update interval, then the performance measures of 
the estimation algorithm are given by the following sample rms values of the estimation 
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errors: 


D ^ ^ 

f^rms — ^ 




- mY 


l-l 


fd 


,rms 


A 1 
~ L \ 


YV^-Ut-)? 


t-i 


(26) 


where L = 4000/TV is the number of update intervals for the entire trajectory. These 
measures are obtained as a function of P/N 0 , where P denotes the received carrier power 
and N 0 is the one-sided power spectral density of the receiver bandpass noise. 

At lower range of (P/No) ratio, the receiver may lose frequency lock, in that the 
frequency errors at times may exceed ± one-half of the sampling frequency f, or ± 
250 Hz. Since the error signals e(Jfc) of (11) are the same for frequency errors of A 
Hz as for the case of A + nf, Hz for any signed integer n, the estimator may make 
frequency estimation errors of nf, Hz. This situation may be referred to as cycle slipping 
in the frequency estimator and is akin to the phenomenon of cycle slipping (phase 
errors equal to multiples of 2tt) in the phase-estimators. If there are one or more cycles 
slipped in frequency, the computed value of / d , rm , would be much larger compared to 
the case when no such cycle slips occur and would be unacceptable. Thus, another 
important parameter for the performance is the probability of maintaining frequency 
lock throughout the trajectory denoted P(lock) or the probability of losing the lock 
P L = l - P(lock). For the purposes of estimating the probability by digital computer 
simulations, 100 simulation runs are made for each value of P/N 0 of interest and an 
estimate of P L is plotted vs the carrier power to noise power spectral density ratio 
(CNR). The sample rms values of (26) are also averaged over all those simulation runs 
for which the frequency lock is maintained. It may well be that for sequences that remain 
under frequency lock, there may be slipping of cycles in the phase estimates. However, 
even under the presence of such cycle slips, the computation made on the basis of (25) 
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provides a good estimate of the pseudo-range as evidenced by the simulations. One 
cycle slip only corresponds to an error of about 0.2 meters in the pseudo-range estimate. 

Figures 6-15 present the simulation results for the DLS algorithm and the composite 
DLS-EKF algorithms presented in the previous sections. The results for the EKF algo- 
rithm operating by itself are available in [14] for comparison. For the simulation results 
a value of A equal to 0.97 has been selected. The initial covariance matrix P(0) for the 
DLS algorithm is selected to be a diagonal matrix, with its diagonal elements equal 
to 2 x io 3 , 2 x io 7 and 2 x io 9 respectively, reflecting the possible uncertainty about the 
parameters. Three different VCO update intervals equal to 5, 10 and 20 sample times 
have been considered. The matrix Q of (19) is also selected to be a diagonal matrix for 
convenience, with its consecutive diagonal elements equal to 4 x io 3 , 2 x io 6 and io 8 . The 
Q matrix represents possible variations in the input signal parameters over an update 
interval and is arrived at from the consideration of a-priori estimate of the maximum 
possible value of the highest order derivative (jerk) present in the input trajectory. For 
the second stage EKF algorithm, the initial covariance matrix £*(0) is selected to be 
also a diagonal matrix but with its elements smaller in value than the corresponding 
elements of the Q matrix, thus effectively resulting in a higher averaging period and 
smaller estimation errors compared to the DLS algorithm. The selected values of diag- 
onal elements of £*(0) matrix are equal to l.o, io 3 and io 6 respectively, in the following 
simulations. From the simulations it appeared to be advantageous, in terms of numer- 
ical stability, to periodically reset the covariance matrix P of the DLS algorithm to its 
initial value. Such a period was selected to be 10 times the VCO update interval. 

Figures 6 and 7 present the simulation results for the performance of DLS algorithm 
while tracking the high dynamic trajectory of Figure 5 in the absence of any data mod- 
ulation. Figure 6 plots the probability of losing the frequency lock P L as a function 
of CNR for two different values of N equal to 10 and 20. As may be inferred from 
the figure, a value of Pl of less than 0.1 is obtained for CNRs above 23.1 dB which 
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is defined to be the threshold point of the algorithm. Figure 6 also plots the average 
number of cycles slipped in the frequency estimation, denoted by W et . In defining such 
a cycle slip, the entire frequency range is divided into disjoint segments of f, Hz with 
the first segment extending from -/,/ 2 to /,/ 2 Hz. Whenever the frequency estimation 
error jumps from one such segment to an adjacent one in either direction, a cycle slip is 
said to occur. Figure 7 plots the rms error in the Doppler estimation as computed from 
(26) and averaged over all convergent sequences. For a CNR between 25 and 30 dB-Hz, 
an rms error of 10-20 Hz is obtained that corresponds to a velocity tracking error of 2-4 
m/s. Figure 7 also plots the average length L e , of a slipped cycle in terms of number 
of samples. The information about J7 C , and T c , is relevant in the case of multistage 
algorithm. Figures 8 and 9 present the results for the probability of losing lock P L and 
the rms estimation error for the DLS algorithm in the presence of data modulation for 
three different values of N equal to 5, 10 and 20. As may be observed from the figures, 
the presence of data modulation increases the threshold by only 0.25-0.5 dB compared 
to the case of no data modulation. The increase in rms frequency estimation error is 
about 10% due to data modulation. 

Figures 10 and 11 present the performance of the composite DLS-EKF algorithm in 
the absence of data modulation. We observe that corresponding to N = 5 , the threshold 
of the algorithm is 22.75 dB-Hz which is slightly smaller than for the DLS algorithm. 
However, the rms estimation errors are significantly smaller than for the single-stage 
DLS algorithm. For the CNR range of 25-30 dB-Hz, the rms error in the Doppler 
estimation lies in the range of 4-15 Hz corresponding to a velocity estimation error 
range of 0.8 to 3 m/ s. The DLS-EKF algorithm also provides the carrier phase estimate. 
The modulo-2x phase-estimation error is plotted in Figure 12, from which it is clear that 
the algorithm is capable of coherent data detection with small probability of error if the 
CNR is higher than 25 dB-Hz. In fact, as shown in Figure 15, the algorithm provides 
good estimates of pseudo-range (related to the absolute phase error via (25)) up to a 
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CNR of about 23 dB at which point the rms error is less than 4 m. The rms pseudo 
range error is less than 1 m for CNRs higher than 25.5 dB-Hz. 

The corresponding results for the performance of DLS-EKF algorithm in the presence 
of data modulation are presented in Figure 13-15. For this case a minimum threshold 
of 23.8 dB is obtained for N = 10 which is about 1 dB higher than for the case of no 
data modulation. In terms of rms estimation errors, for a CNR range of about 25-30 
dB-Hz, the rms frequency estimation error lies in a range of 8-20 Hz corresponding to 
a velocity error of about 1.5 to 4 m/s. For this case, as is apparent from Figure 15, the 
pseudo-range estimation errors are also higher and for a CNR range of 25-30 dB-Hz lie 
in a range of 0.3-6 m. Notice, however, that no sharp threshold is observed in either 
the frequency or phase estimation errors over the entire range of CNR between 22-30 
dB-Hz considered in the simulations. 
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7. COMPARISON WITH PREVIOUS TECHNIQUES 

For the case of no data modulation, we compare the performance of the proposed al- 
gorithm with some of the techniques analyzed in [14] in terms of their performance when 
tracking exactly the same high dynamic trajectory. Compared to a more computation- 
intensive maximum likelihood estimate (MLE) of [14], we observe that the DLS-EKF 
algorithm requires about 0.25 dB smaller CNR than MLE in terms of threshold. In 
terms of rms frequency estimation errors, the MLE achieves an rms error between 8 Hz 
to 35 Hz at a CNR of 23 dB-Hz depending upon the estimation delay ranging between 
30-80 samples (higher delay provides smaller error). The DLS-EKF algorithm provides 
an error of 35 Hz for a delay of 5 samples at a CNR of 23 dB-Hz. The MLE algorithm 
does not provide any phase estimate. Compared to single stage EKF algorithm, we 
observe that the DLS-EKF algorithm is about 1. 5-2.0 dB better in terms of threshold 
depending upon the value of exponential data weighting coefficient and the filter order 
used in the simulations of [14]. In terms of RMS errors, the performance is similar to 
that of third order EKF alone. Notice, however, that direct comparison with the results 
of [14] may be somewhat misleading. This is so because while the DLS-EKF algorithm 
includes all of the sequences in the computation of rms error above 25.5 dB-Hz, EKF 
rejects about 5% of the worst sequences, as the probability of losing lock is about .05 at 
CNR of 25.5 dB-Hz. For the case of cross product AFC loop of [14], the threshold lies in 
a range of 25-28 dB-Hz depending upon the loop parameters. Thus, the DLS-EKF algo- 
rithm is superior by 2-5 dB-Hz compared to AFC loop. AFC loop provides a minimum 
rms frequency error of 25 Hz at a CNR of 28 dB-Hz compared to a minimum of 5 Hz 
achieved for the DLS-EKF algorithm for the same CNR. Notice that in an AFC loop, the 
parameters achieving a relatively low estimation error are different than those yielding 
low thresholds and thus, a range of loop parameters must be considered for proper 
comparison. In terms of rms phase error the performance of the DLS-EKF algorithm is 
similar to EKF alone. In terms of computations the DLS-EKF algorithm of the publi- 
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cation requires about the same number of computations as for a fourth-order EKF but 
about twice as many computations as a third-order EKF. The number of computations 
are at least an order of magnitude smaller than the MLE. 

For the case when the data modulation is present, we compare the performance of 
the DLS-EKF algorithm with the MLE algorithm of [7], where a somewhat less severe 
GPS trajectory is analyzed. The results of [7] show a marked threshold of about 30 
dB-Hz in terms of CNR compared to 24 dB-Hz for the proposed algorithm. Thus, 
the proposed algorithm results in a 6-dB reduction in terms of threshold compared 
to previous schemes of the literature. Above threshold both the DLS-EKF and MLE 
algorithms have similar values for rms range estimation and range rate estimation errors 
(minor differences result due to differences in the trajectories). In terms of computations 
both of the algorithms are comparable. In terms of threshold on CNR, the DLS algorithm 
is very close to the composite DLS-EKF algorithm. However, in terms of rms frequency 
estimation errors, it has significantly higher estimation errors. In those cases where 
higher estimation errors are acceptable, one may apply the DLS algorithm by itself, as 
it requires only one-half of the computations required by the DLS-EKF algorithm. 



8. CONCLUSIONS 


This publication has presented a novel multistage estimation scheme for the effi- 
cient estimation of the phase and frequency of a very high dynamic signal, which may 
possibly be phase modulated by unknown binary data and is received under relatively 
low carrier-to-noise power ratio conditions. The proposed scheme is of very general 
nature and has much wider scope than the applications in this publication. For a very 
important application considered in this publication, the algorithm has been special- 
ized to have just two stages. The first stage of the estimation scheme is a least-squares 
algorithm operating upon the differential signal model while the second stage is an 
extended Kalman filter of third order. For the very high dynamic GPS trajectories, 
the proposed algorithm has been shown to significantly outperform the previous algo- 
rithms of the literature in one or more aspects, including threshold on CNR, estimation 
errors, availability of phase estimates and thus the estimate of pseudo-range, compu- 
tational complexity, and flexibility. For the case of no data modulation, the proposed 
scheme has a threshold that is slightly lower than the more computation-intensive im- 
plementation of MLE algorithm. When compared to just the EKF operating by itself, 
the proposed DLS-EKF scheme provides from about 1.5- to 2-dB reduction in threshold. 
In comparison to more conventional schemes, such as AFC loops, the performance is 
even better. 

For the case when an unknown data modulation is present, the algorithm provides 
an improvement of 6 dB in terms of threshold on CNR in comparison to the MLE scheme 
of [7] specifically proposed for such applications. In addition to phase and frequency 
estimates, the algorithm can provide estimates of frequency derivative as well, although 
not presented in this publication. The scheme being of a very general nature, it may 
be possible to reduce the threshold even further by using a higher dimension for the 
state vector related to the higher number of terms in the Taylor series expansion in 
arriving at the signal model for the first stage DLS algorithm. Further improvements 
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are possible by the application of more optimum sampling techniques as proposed in 
[11]. The performance may also be improved both in terms of the threshold and the 
rms estimation errors by increasing the number of stages to three or more. 
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Figure 1. Signal Parameters Estimation by DLS Algorithm 
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Figure 2. LS Algorithms (a) LS Algorithm With Model Noise Colored (b) LS Algorithm With 
Prewhitened Noise - Conceptual (c) An Approximate and Realizable Equivalent of (b) 
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Figure 3. Single-Stage Estimators (a) A Generalized Single-Stage Estimator 

(b) An Equivalent Form of Single-Stage Estimator 
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Figure 4. A Multistage Estimator for the Process 0i(t) 
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s 7. RMS Frequency Estimation Error vs CNR; DLS Algorithm in the Absence of 
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Figure 9. RMS Frequency Estimation Error vs CNR; DLS Algorithm in the Presence of 
Data Modulation 
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Figure 10. Probability of Losing Frequency Lock vs CNR; DLS-EKF Algorithm in the Absence of 
Data Modulation 


37 



RMS FREQUENCY ESTIMATION ERROR (Hz) 



CNR (dB-Hz) 


Figure 11. RMS Frequency Estimation Error vs CNR; DLS-EKF Algorithm in the Absence of 
Data Modulation 
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Figure 12. Modulo 2 ji RMS Phase Error vs CNR; DLS-EKF Algorithm in the Absence of 
Data Modulation 
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Figure 13. Probability of Losing Frequency Lock vs CNR; DLS-EKF Algorithm in the Presence of 
Data Modulation 
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Figure 14. RMS Frequency Estimation Error vs CNR; DLS-EKF Algorithm in the Presence of 
Data Modulation 
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Figure 15. RMS Pseudo-Range Estimation Error vs CNR for DLS-EKF Algorithm; (With and Without 
Data Modulation) 
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